rm(list=ls())
library(foreign)
library(locfit)
library(Hmisc)
setwd("/home/users/innoshu.AD3/drdfinal")

setwd("/Users/innoshuadmin/Dropbox/1RD/revision/drdfinal")
source("./functions/DRDfunctions.R")
## load in robust bandwidth selection functions from CCT(2014)## 
source("./functions/rdbwselect.R")
source("./functions/functions.R")

source("./APP3/run_ab.R")


dat <- read.dta(file="./APP3/app3.dta")
dat$def=dat$lncn*0
dat$def[dat$anno==1993]=log(0.91324)
dat$def[dat$anno==1998]=log(1.07596)
dat$def[dat$anno==2000]=log(1.12095)
dat$def[dat$anno==2002]=log(1.178967)
dat$def[dat$anno==2002]=log(1.232134)
dat$lncn=dat$lncn-dat$def
dat$lncn_n=dat$lncn_n-dat$def
dat$lncn_rtn=dat$lncn_rtn-dat$def
dat$lnfood=dat$lnfood-dat$def
dat$lnfood_n=dat$lnfood_n-dat$def
dat$lnfood_rtn=dat$lnfood_rtn-dat$def

newdata<- dat[c("lncn","lncn_n", "lncn_rtn","lncn_oecd", "lncn_rtoecd", "lnfood","lnfood_n", "lnfood_rtn","lnfood_oecd", "lnfood_rtoecd","esse_m","esse_f","pop_category","qu_m","famtype","strata","Hs")]
newdata$qu_m=(newdata$qu_m=="retired")*1
newdata$famtype=(newdata$famtype=="couple")*1

H=0
uniquestrata=unique(newdata$strata)
NStrata=length(uniquestrata)
for (vstrata in 1:NStrata){
H =H+mean(newdata$Hs[newdata$strata == uniquestrata[vstrata]])
}

RepeatTimes<-runif(202)*0
for (vstrata in 1:NStrata){
RepeatTimes[uniquestrata[vstrata]] = length(newdata$strata[newdata$strata == uniquestrata[vstrata]])
}
newdata$Ns = RepeatTimes[newdata$strata]
N=nrow(newdata)

newdata <-subset(newdata,is.na(lncn)==0&is.na(lnfood)==0&esse_m<=20& esse_m>=-30)
esse_f_nan<-is.na(newdata$esse_f)

Xeval <- seq(-30, 20, length=21)  
ctff<-0
################################################################

print("Table 7")
for (numy in 1:6){
	if (numy==1) {yname="lncn"} else
	if (numy==2) {yname="lnfood"} else
	if (numy==3) {yname="lncn_rtoecd"} else
	if (numy==4) {yname="lnfood_rtoecd"} else
	if (numy==5) {yname="lncn_oecd"} else
	if (numy==6) {yname="lnfood_oecd"} 
numd=1	
results<-run_ab(numd=numd,yname=yname)
print(yname)
print("Proportion of Compliers (unadjusted-- sampling weight adjusted--#obs)")
print(round(results$pc,3))
print("Mean RD (unadjusted-- sampling weight adjusted--#obs)")
print(round(results$mrd,3))
}


for (table in 1:4){
if (table>=3){
print(paste("Robustness Check"))
} else{	
print("Table 8")
}
for (numy in 1:2){
	if (numy==1) {yname="lncn_rtoecd"} else
	if (numy==2) {yname="lnfood_rtoecd"} 
numd=4		
results<-run_ab(numd=numd,yname=yname,table=table)
print(yname)
print("Proportion of Compliers (unadjusted-- sampling weight adjusted--#obs)")
print(round(results$pc,3))
print("Mean RD (sampling weight adjusted--#obs)")
print(round(results$mrd[,4:7],3))
print("Distributional RD (sampling weight adjusted--#obs)")
print(round(results$drd[,4:7],3))
}
}





